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ABSTRACT 

HERO (Hybrid Evaluator for Radiative Objects) is a 3D general relativistic radiative 
transfer code which has been tailored to the problem of analyzing radiation from 
simulations of relativistic accretion discs around black holes. HERO is designed to be 
used as a postprocessor. Given some fixed fluid structure for the disc (i.e. density and 
velocity as a function of position from a hydrodynamics or magnetohydrodynamics 
simulation), the code obtains a self-consistent solution for the radiation field and for 
the gas temperatures using the condition of radiative equilibrium. The novel aspect 
of HERO is that it combines two techniques: 1) a short characteristics (SC) solver 
that quickly converges to a self consistent disc temperature and radiation field, with 
2) a long characteristics (LC) solver that provides a more accurate solution for the 
radiation near the photosphere and in the optically thin regions. By combining these 
two techniques, we gain both the computational speed of SC and the high accuracy 
of LC. We present tests of HERO on a range of ID, 2D and 3D problems in flat space 
and show that the results agree well with both analytical and benchmark solutions. 
We also test the ability of the code to handle relativistic problems in curved space. 
Einally, we discuss the important topic of ray-defects, a major limitation of the SC 
method, and describe our strategy for minimizing the induced error. 

Key words: methods: numerical - radiative transfer - accretion, accretion discs - 
black hole physics 


1 INTRODUCTION 

Radiative transport (RT) plays a crucial role in astrophysics. 
It governs our primary source of information about the cos¬ 
mos, viz., the luminosities and spectra of cosmic sources^. 
In addition, radiation acts as a channel for energy trans¬ 
port, shaping the dynamics and impacting the evolution of 
astrophysical systems on all scales: stellar evolution, stel¬ 
lar and galactic winds, super-Eddington accretion, epoch of 
reionization, etc. 

Radiation is quite challenging to model since it behaves 
in highly nonlinear and nonlocal ways. The radiation field 
is determined by a six-dimensional (6D) system of coupled 
equations which link spatial positions (3D), ray directions 
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(2D) and frequencies (ID), all as a function of a seventh 
coordinate, time (in the general time-dependent problem). 
The high dimensionality of the solution vector, combined 
with the complex integro-differential nonlocal nature of the 
problem, results in an extremely taxing computational chal¬ 
lenge (both in terms of memory and computational speed). 

Because of the intrinsic complexity of the problem, the 
earliest radiative solvers made use of fairly restrictive ap¬ 
proximations. For example, enforcing spatial symmetries al¬ 
lows one to greatly reduce the dimensionality of the problem 
and led to the first generation of 1D/2D radiative codes (e.g. 
Feautrier method that exploits symmetric and antisymmet¬ 
ric radiation moments, see Mihalas et al. 1978). Several of 
the multidimensional RT ideas have also found application 
in the field of neutrino transport in core collapse supernovae 
(Burrows et al. 2000; Rampp & Janka 2002; Liebendorfer et 
al. 2004; Livne et al. 2004; Hanke et al. 2013). In this case, 
the infalling matter becomes optically thick to the neutrino 
flux, which propagates analagously to radiative transport in 
optically thick media. 

3D radiative problems are extremely computationally 
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taxing, which motivates the need for highly efficient tech¬ 
niques. Recent applications that make use of 3D RT in¬ 
clude models of young stellar objects (Wolf, Fischer Sz Pan 
1998), protostellar to protoplanetary discs (Indebetouw et 
ah 2006; Niccolini & Alcolea 2006), reflection nebulae (Witt 
& Gordon 1996), molecular clouds (Steinacker et ah 2005; 
Pelkonen, Juvela & Padoan 2009), spiral galaxies (Bianchi 
2008; Schechtman-Rook, Bershady & Wood 2012), interact¬ 
ing and starburst galaxies (Chakrabarti et ah 2007; Hayward 
et ah 2011), AGNs (Schartmann et ah 2008; Stalevski et ah 

2012) , and cosmological reionization simulations (Ferland et 
ah 1998; Abel & Wandelt 2002; Finlator et ah 2009). These 
applications are primarily concerned with point sources il¬ 
luminating optically thick media. 

In more complex optically thick flows where extended 
diffuse emission is important, the typical approach is to sim¬ 
plify the directional structure of the radiation field. One idea 
is to decompose the radiation field into its moments (the first 
three moments being the radiation energy density, radiation 
flux, and radiation pressure) and to evolve the radiation field 
locally as a fluid using some closure relation on the moments. 
This approach is popular due to its local nature and fast 
speed, and has been implemented in many hydrodynamic 
codes (Turner & Stone 2001; Bruenn et al. 2006; Hayes et 
al. 2006; Gonzales, Audit, & Huynh 2007; Krumholz et al. 
2007; Gittings et al. 2008; Ohsuga et al. 2009; Swesty & 
Myra 2009; Gommercon et al. 2011; Zhang et al. 2011; Kolb 
et al. 2013). Flux limited diffusion (FLD, Levermore & Pom- 
raning 1981), which is based on the first two moments, is 
the simplest and perhaps most popular moment-based ra¬ 
diative method. Ml closure is a generalization of the FLD 
method which uses the first three moments. It has recently 
gained traction as a fast solver for diffuse emission (Lever- 
more 1984; Dubroca & Feugeas 1999; Stone et al. 1992; 
Gonzales, Audit, & Huynh 2007; Sadowski et al. 2013). 

The field of protoplanetary disc dust modelling is one 
area that has been active in developing RT techniques. In 
this domain, the monte-carlo approach has been the domi¬ 
nant paradigm owing to its ease of implementation and abil¬ 
ity to handle anisotropic scattering kernels (Lopez, Mekar- 
nia, & Lefevre 1995; Niccolini, Woitke, & Lopez 2003; Wolf, 
Henning & Stecklum 1999; Bjorkman & Wood 2001; Pinte, 
Duchene, & Bastien 2006). Monte-carlo RT has also been ap¬ 
plied to the problem accretion discs, specifically to handle 
the problem of modeling Gompton scattering in disc coro¬ 
nas (Dolence et al. 2009; Schnittman & Krolik 2013; Ghosh 

2013) as well as dust scattering of emission from the Galac¬ 
tic centre (Odaka et al. 2011). The primary downside of 
this technique is the computational expense and noise asso¬ 
ciated with photon statistics. The method is also limited to 
modest optical depths. 

The limitations of the various methods mentioned above 
led to the development of more deterministic discretized 
finite-differencing methods (Stenholm, Stoerzer, & Wehrse 
1991; Steinacker, Bacmann, & Henning 2002) and dis¬ 
crete ordinates characteristic methods (Dullemond & Tur- 
olla 2000; Steinacker, Bacmann, & Henning 2002; Hayes & 
Norman 2003; Vogler et al. 1982; Heinemann et al. 2006; 
Woitke, Kamp & Thi 2009; Hayek et al. 2010; Davis et 
al. 2012; Jiang et al. 2014), which have a number of widely 
discussed advantages. However, all of these methods suffer 
from the phenomenon of ” ray-defects” (see appendix A for 


a detailed discussion). Another twist to the discretization 
strategy involves expanding and evolving the radiation field 
as a combination of spherical harmonics (Szu-cheng & Kuo- 
Nan 1982; Evans 1997; McGlarren, Holloway, Sz Brunner 
2008). Spherical harmonics respect rotational symmetry and 
therefore do not suffer from the linear ray defect patterns 
that plague other discrete ordinate methods. The tradeoff is 
that spherical methods typically suffer from artificial side- 
lobe patterns due to the finite order cutoff used in taking 
the series expansion of the radiation field. 

Raytracing methods are also used for predicting the ob¬ 
served flux from astrophysical objects. These codes typically 
do not solve for the global radiation field within an object; 
instead the focus is on calculating the apparent intensities 
that reach a distant observer. In these codes, the typical 
assumption is to ignore scattering (i.e. to ignore nonlocal 
coupling of the radiation field) since this allows one to im¬ 
mediately compute the evolution of ray intensity by simply 
integrating the local emissivities along a photon geodesic 
(e.g. Gunningham & Bardeen 1973; Ozel & Di Matteo 2001; 
Huang et al. 2007; Dexter & Agol 2009; Shcherbakov & Lei 
2011; Vincent et al. 2011; Ghan et al. 2013; Bohn et al. 2014). 

Despite the multitude of radiative solvers available, 
none are currently able to solve the problem of optically 
thick emission from accretion discs around black holes 
(BHs). The main difficulty here is that a general relativis¬ 
tic 3D framework is needed to properly account for both 
light bending effects and doppler/gravitational redshifts. 
The work reported here is a first attempt at tackling the 
radiation problem in full glory around black holes. The code 
we describe here called HERO is a postprocessor - given 
some fixed gas structure (density, velocity, energy injection 
rate), as determined by a separate BH accretion disc simu¬ 
lation, our code solves for the radiation field along with its 
self-eonsistent temperature solution. We explain here how 
our code works and show verification tests to gauge its per¬ 
formance under various conditions. 

The organization of the paper is as follows. In §2, we 
first describe the method by which HERO operates. We 
give a brief overview of the radiative transfer problem that 
HERO solves, and how it generalizes in curved space. We 
also detail how the two radiative solvers (short/long char¬ 
acteristics) operate, and the numerical methods involved 
(i.e. acceleration schemes, interpolation, discretization strat¬ 
egy). Next, in §3 we check the code by comparing to ana¬ 
lytic and benchmark results for ID, 2D, and 3D test con¬ 
figurations. Einally, in appendix §A we end with a discus¬ 
sion of ray-defects, a systematic problem that plagues short- 
characteristic radiative solvers. 


2 RADIATIVE SOLVER 

Our radiative code consists of three primary components: 

(i) A short characteristic solver for quickly obtaining 
an approximate solution to the radiation field 

(ii) A long characteristic solver for more accurate 
modeling of the radiation field 

(hi) An optional raytracer for generating mock observa¬ 
tions from some distant observing plane 
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We employ a hybrid approach to solve the radiation 
field, hence the name of the code: HERO (Hybrid Evaluator 
for Radiative Objects). The code is hybrid in the sense that 
it uses both (i) “short characteristics” (SC, Mihalas et ah 
1978) and (ii)“long characteristics” (LC, Eeautrier 1964). 
Here “short” and “long” refer to the length of the light rays 
that are considered in each iteration of the solver. The short 
characteristics method is only concerned with propagation 
of radiation from a given cell to its immediate neighbours, 
whereas the long characteristic method traces rays all the 
way to the edge of the computational grid. 

The motivation for developing a hybrid approach is 
to allow for the accurate modelling of radiation in opti¬ 
cally thin regions via long characteristics, while retaining 
the computational speed offered by short characteristics in 
optically thick regions (see §2.1 & §2.3 for a detailed com¬ 
parison of the methods). Eor any given problem, we first 
apply short characteristics to quickly solve for the local ra¬ 
diation/temperature. We then feed this solution as the ini¬ 
tial guess for a more detailed long characteristics calculation. 
Using such a hybrid approach allows us to combine the “best 
of both worlds.” 

In addition to solving for the radiation field, in both SC 
and LC, we also solve for a self-consistent gas temperature. 
The general idea is to loop back and forth between solv¬ 
ing for the radiation field given a fixed temperature struc¬ 
ture, and solving for the equilibrium temperature distribu¬ 
tion given a fixed radiation field. This procedure is iterated 
until convergence, leaving us with a self-consistent solution 
for both the radiation and the temperature. 

Our goal in developing HERO is to model/investigate 
the observed properties of relativistic accretion discs around 
black holes, and to compare the radiative properties of sim¬ 
ulated discs with data collected by earthbound telescopes. 
Therefore, after completing the SC and LC steps described 
above, the final step is to input the radiative and temper¬ 
ature structure as computed from LC and to generate syn¬ 
thetic observations of the disc as seen by a distant observer. 
Eor this stage we use the standard raytracing approach (see 
§2.5 for details). 

Due to the high temperatures present in accretion discs, 
Comptonization plays a crucial role in determining the ob¬ 
served radiation from relativistic disc photospheres. Because 
of the complex nature of the Compton scattering kernel, 
we defer discussion of Comptonization to a follow-up paper, 
which will focus exclusively on explaining and testing our 
relativistic scattering module. 


2.1 Short Characteristics 

Short characteristics is a popular algorithm for tackling mul¬ 
tidimensional radiative transfer problems (Mihalas et al. 
1978; Olson & Kunasz 1987; Kunasz & Auer 1988). It is 
a local method and hence very fast. The basic idea is to 
solve for the radiation field at a given point (reference cell) 
using only the information provided by neighbouring cells 
(see schematic in Eigure 1). 
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Figure 1. Schematic of the short characteristics method. A 
curved ray (null geodesic) is shot “upstream” to determine the 
intensity at F. This ray is terminated at the neighbouring cell 
boundaries, and the intensity Iq at the boundary is computed by 
interpolation of neighbouring grid points (in this case points H 
and I). The source function S(t') in Eq. 22 is also obtained by 
interpolation on neighbouring points (points E, F, H, I), and the 
intensity at F is thereby calculated. The procedure is repeated 
for a number of rays in different directions to obtain an estimate 
of the radiation field at F. 


2.1.1 Ordinary Radiative Transfer Equation 

To set the stage, we first discuss the standard case of flat 
space, where the radiative transfer equation takes the form: 

+ cru)Iu + ju + CTu J 0i.(U,U')A(Q')dU'. (1) 

Here E is the intensity of a ray travelling along path s di¬ 
rected towards U; hii, and are the absorption and scatter¬ 
ing coefficients; (j)jy is the scattering shape function, normal¬ 
ized such that f (t)jy(Tl,Tl')d^' = 1; and jjy is the emission 
coefficient. Typically, Eq.l is simplified to the form 


Q 

, — J-jy \ 

di~i/ 

where the the optical depth Tu is given by 
dru = {k,u + crjy)ds, 
and the source function Si^ is defined as 

cr _ f (t>u{0,Q')E{0')dO' 

Ou — I . 

Kii/ — 0^1/ 

In the case of isotropic scattering, = l/47r, and the scat¬ 
tering term simplifies to 


( 2 ) 

(3) 

(4) 




(5) 


where Ju is the angle-averaged mean intensity. We assume 
isotropic scattering in all of the work reported here. We fur¬ 
ther assume a thermal source, for which the emissivity is 
given by ji, = KuBu{T), where Bjy{T) is the Planck function 
corresponding to the local temperature T of the medium. 
With these two simplifications, the source function can be 
compactly rewritten as: 


Su — ejyBi, + (1 ej^)Jjy, 


(6) 
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where is the photon interaction destruction probability, 
given by the ratio of the absorption and total opacity: 


Kjy 

Cjy = ■ . 

kljy + CTjy 


(7) 


To avoid confusion with summation indices, in the remainder 
of this paper we will no longer use a subscript iz to note the 
frequency dependence of radiative quantities. Unless oth¬ 
erwise stated, this will apply to all subsequent intensities, 
opacities, source functions, and emissivities. 


2.1.2 Covariant Formulation of the Radiative Transfer 
Equation 

Because HERO has been developed for relativistic problems, 
we use a generalization of the radiative transfer equation 
that accounts for relativistic effects such as light bending 
and doppler boosting. The generalization of the radiative 
transfer equation to curved space has been discussed in the 
literature (Mihalas & Mihalas 1984; Lindquist 1966), but we 
give a brief primer here for completeness. Starting from the 
time dependent version of Eq. 1: 

if + (*.V)/ = G, (8) 

where G on the righthand side contains all scattering and ab¬ 
sorption source terms that interact with the radiation field, 
and n is a normalized 3-vector denoting the propagation 
direction of photons. We can rewrite Eq. 8 in a more re¬ 
vealing form by introducing the direction 4-vector whose 
orthonormal form is (note that k - k — C)\ 

r = (l,n). (9) 

This lets us rewrite the RT equation (setting c— 1) as: 

= G = -(n + a)I +j-\-(jJ. (10) 

By considering the Lorentz transformation properties of 
the various quantities, Eq.lO can be recast in an invariant 
form (see Mihalas & Mihalas 1984, §7.1): 

where the photon 4-momentum = hizk^ takes on the role 
of the propagation 4-vector. Instead of the usual definitions 
of intensity I and source term G, this equation considers 
their corresponding relativistically invariant versions, 

I=Ilv^, g = hGlv^, ( 12 ) 

where the scaling of / is determined by photon number con¬ 
servation in phase space. The (hiz)~‘^ scaling for G arises 
after transforming Eq.lO to use X and allows us to construct 
the invariant scalings for the absorption and scattering co¬ 
efficients as: 


II 

, 5i/ — hizo'i/^ 

(13) 

and for the emissivity and 

= hj^/v 

source function, 

^ S^S/iz\ 

(14) 


In terms of these new quantities, the invariant source Q be¬ 
comes 


g^-{lFB){X-S). (15) 


Eor convenience, since we usually work with comoving quan¬ 
tities, /^c, cTc, aS'c, we can instead write Q as: 




— hlZc(kic + CTc) 



= + CTc) 



(16) 


where is the four-velocity of the fluid. Eq. 11 is in mani¬ 
festly covariant form and the generalization to curved space 
is simply accounted for by replacing the directional deriva¬ 
tive dfdx^ with the covariant derivative: 

Dx^ dx^ ^ ^ 

evaluated along a null geodesic (Lindquist 1966). The gen¬ 
eral covariant form of the RT equation then takes the form: 

= ^ (18) 

We can further simplify the above expression by using the 

photon geodesic equation, 

^ (19) 

where A is an affine parameter defined by = dx^ f d\. 
Using Eq. 19 to eliminate the Christoffel symbols in Eq. 17, 
the RT equation becomes: 


e, dX dp^ dX 

^ dx^ dX dp^ 




( 20 ) 


Einally, using the chain rule for partial differentiation, we 
obtain the very simple equation 


dX 

dX 




( 21 ) 


This is the form of the covariant RT equation that we 
employ in HERO. The only point that requires some care is 
choosing the correct frequency when evaluating the source 
term Q. Eor instance, if we wish to evaluate X at frequency 
iz in the lab frame, but we compute Q in the fluid comoving 
frame (the natural frame since the opacity is defined there), 
then iZc in Eq. 16 must be chosen such that it corresponds 
to the lab frame iz used to label X. 


2.2 Implementation of Short Characteristics 

Given a ray in the reference cell E, we compute the up¬ 
stream location where it intersects the boundary of that 
cell (see Eig. 1) and use this location to determine the rela¬ 
tive weights contributed by each of the four nearest neigh¬ 
bours on that boundary. In flat space, it is possible to write 
down an analytical formula for the intersection point (e.g. 
Dullemond & Turolla 2000), but there is no generalization 
in terms of simple functions for curved space. HERO does 
the calculation numerically. 

Although this is computationally somewhat expensive, 
it needs to be done only once as part of the initialization 
steps. The information from this initial computation is saved 
and used repeatedly during the SC iterations. Given the in¬ 
tersection of a ray with the nearest neighbour cell boundary, 
we compute the intensity of this ray in the reference cell E 
using the radiative transfer equation (Eq. 21): 

/f =/oexp(-To) + f S{r')e:^p{-T')dr', (22) 

Jo 
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where Iq is the incoming intensity at the boundary, and tq 
is the total optical depth of the photon path leading from 
the boundary to point F (computed using the average of the 
central and boundary opacities). The intensity Iq is com¬ 
puted by interpolating from the four boundary points. Our 
current implementation of short characteristics uses linear 
interpolations at various stages (this is easily improved in 
the future). In this spirit, the integral evaluation in Eq. 22 
is computed as 

S{r') e:>cp(—r')dr' 

_ e ^{So — Sr{lr))So{r — 1)Sr ^^^ 3 ^ 

which is the analytic solution when the source function S 
varies linearly with r along the photon trajectory. Higher or¬ 
der interpolation schemes have been explored in other codes 
(e.g. Kunasz & Auer 1988); however, special care is needed 
to ensure non-negativity of the source function. 

Our implementation of SC in HERO retains some mem¬ 
ory of the radiative quantities from the past iteration. This 
helps with stability. Thus, in computing the intensity at it¬ 
eration n, we use the linear combination 

7" = (1 - m)I^s^ + mr-\ (24) 

where is the result from the RT formal solution (i.e. 

evolving the radiation field through SC). Typically, we set 
m = 0.5. HERO with SC converges within a reasonable 
number of iterations ~ 10 — 1000 s, depending on the degree 
of scattering/coupling in the problem. 

The above discussion assumed that the gas tempera¬ 
tures are given. More often, one does not know the temper¬ 
ature but has to solve for it based on boundary conditions 
and internal heat sources. In this case, between every two 
of the above iterations for the intensity, the code carries out 
a round of iterations (typically ~ 10 s) to solve for the tem¬ 
perature using the condition of radiative equilibrium, i.e., 
the requirement that the heating and cooling rates of the 
gas should balance. Radiative equilibrium requires at each 
position X 

Qcool(^) — Qheat(^) T Qinject(^)5 

where 

Qcooi(^c) = J hij,(x)Bjy[T{x)]diy , 

Qheat(^) — l^i/{x^ Ji/{x^dlO^ (26) 

and Qinject(a^) refers to any additional injection of energy 
into the fluid, e.g., through viscous dissipation in an accre¬ 
tion disc. 



2.3 Long Characteristics 

Similar to the short characteristics method, our long char¬ 
acteristics solver obtains the radiation field by evolving ray 
intensities according to the radiative transfer equation (c.f. 
Eq. 21 ). Eigure 2 shows a schematic of our approach to LC 
- rays are shot upstream, and we evaluate the intensity at 
the observing cell according to 

/o = /6exp(-T6) + f S{r')e:^p{-r')dr'. (27) 

Jo 
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Figure 2. Schematic of long characteristics RT solver. We cal¬ 
culate the local mean intensity at an observing cell (large blue 
circle) by integrating the source function along photon geodesics 
(light blue path). The RT integral is evaluated by direct summa¬ 
tion of contributions from linear piece-wise segments (small blue 
dots). 


The above expression is identical to our SC calculation 
(Eq.22), except that A now represents the intensity distri¬ 
bution at the grid boundaries and the ray integral continues 
past the immediate neighbouring cells to traverse the entire 
spatial domain. The integral in Eq. 27 is evaluated in dis¬ 
crete steps (r r At) , where each source-function piece 
is evaluated using the same linear interpolation scheme as 
used in SC (Eq. 23). 

One major difference between our LC and SC solvers 
is the choice of angular grid. In SC we used a fixed angu¬ 
lar resolution (Na — 80, see §2.7.3) for all cells. This fairly 
low resolution choice is adequate (though only barely) in 
SC since one only needs to resolve the neighbouring 27 grid 
cells. In LC, since the rays traverse the entire spatial do¬ 
main, much higher angular resolution is needed to resolve 
the emission from distant source cells. This becomes partic¬ 
ularly problematic with spherical coordinates where the grid 
spans a few decades in radius. Cells located at larger radii 
have a difficult time “seeing” the inner regions unless one 
takes care in choosing the ray directions. 

Our solution is to dynamically set the angular resolution 
used by the LC solver such that there is sufficient coverage 
to resolve all the cells within the domain. Consider, as a 
typical example, a 3D grid in spherical coordinates which 
extends from an inner radius rin to an outer radius Tout 
Tin. The radial dynamic range may be up to three decades, 
as in some of the examples described later. Let us suppose 
that our spherical spatial grid uses uniform angular spacing 
in with Ne cells, and logarithmic spacing in radius with 
Nr < Ne cells per decade. Consider a cell at some radius r. 
Eor this cell, the apparent density of other cells as viewed 
in its local ”sky” is highly anisotropic. Eor rays traveling 
from the direction of the coordinate centre, this cell needs 
to sample the local radiation field with an angular resolution 
of order AOres ~ Hn/'riV^ in order to adequately sample the 
contribution of all the cells at the centre. In the opposite 
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direction, however, an angular resolution of order ~ I/Nq is 
sufficient. 

The above requirement can be achieved by designing a 
suitable angular tiling of the sky. One simple prescription is 
to vary the angular resolution A^res as a function of direction 
X relative to the local radius vector as follows: 

where / controls the factor by which we oversample to re¬ 
duce noise. Usually / = 1 is adequate, however even with 
this relatively low resolution, the number of tiles varies from 
~ Nq for the innermost radial grid cells to several orders of 
magnitude larger for outer cells. The increase in resolution 
at the outer cells is needed in order to resolve the interior. 
The average number of rays per cell is thus a factor of 100 
or more larger than the Na — 80 rays that we consider with 
short characteristics. In addition, each ray in LC traverses 
the entire grid, compared to just one cell with SC. For both 
reasons, one iteration of LC often takes 10^ times more com¬ 
putation time than one iteration of SC! On the other hand, 
one iteration of LC moves the solution much closer to con¬ 
vergence since information is propagated across the entire 
grid. 

The main weakness of LC is its inability to handle op¬ 
tically thick regions in scattering-dominated problems. This 
is because of the huge number of iterations needed to model 
properly the diffusion process. Since LC cannot be feasibly 
run for more than a few dozen iterations, it cannot be used 
to update radiative information deep inside optically thick 
objects. Luckily, it is precisely in this optically thick regime 
where short characteristics is free from the ray-defect prob¬ 
lem and can be trusted to produce unbiased results. This 
motivates our hybrid scheme: first we run many iterations of 
short characteristics to resolve the complete diffusive process 
and to correctly obtain the radiation field in the optically 
thick parts of the problem; second we run a few passes of 
LC to ensure that the optically thin regions of the radiation 
field are handled accurately and are devoid of ray-defects. 


2.4 Acceleration Schemes 

The radiative quantities S and J are inextricably linked. For 
future reference, we quantify the dependence of J on S' by 
means of a A operator, defined such that 

= (29) 

where the sub/superscripts denote grid cells. Each iteration 
of the radiative solver^ acts as A, describing how J is related 
to S. However, since S itself has a scattering term, it needs 
to be updated along with J. Writing this explicitly (c.f. Eq. 
6 ), 

^U+I = (1 _ e)J’^ + eB 

= (1 - e)A[S"*] + eB. (30) 

By switching between solving for J as a function of S (Eq. 
29), and for S' as a function of J (Eq. 30), we seek to find a 


^ In both SC and LC, the mean intensity term J is evaluated by 
averaging I via Eq.5 


stationary self-consistent solution for all our radiative quan¬ 
tities. Unfortunately, simply applying these two transfor¬ 
mations one after the other and iterating has well known 
convergence issues whenever scattering strongly dominates 
(say with e < 10“^). 

To avoid slow convergence, acceleration schemes are of¬ 
ten employed. Accelerated lambda iteration (ALI) is one 
such technique. To understand the idea behind ALI, rewrite 
Eq. 30 so as to isolate S\ 

5= [l-(l-€)A]“^[eB]. (31) 


Eormally, one can obtain the complete solution by evaluat¬ 
ing the inverse matrix on the right hand side of Eq.3I and 
calculating S from B. Then evaluating Eq. 29 gives J. How¬ 
ever, in 3D, A is an enormous matrix and it is impractical 
to carry out full matrix inversion. Instead, ALI considers 
an approximate lambda operator, typically chosen to be the 
diagonal component Ku (i.e. the contribution to the local 
J from solely the local S - usually the largest component). 
Denoting the shift in our old formal solution without accel¬ 
eration as solving the system taking 

into account A yields (see Hubeny 2003 for details): 


= 


AS 


'FS 


1 - (1 - ei)Ai 


(32) 


To obtain the diagonal of the A operator, we consider: 


A« = ^, (33) 

where Ju represents the local contribution to J from the 
source function within the reference grid cell. This self¬ 
illumination contribution to J is the following sum over the 
NA ray angles considered: 


1 

J _ A ^ ^ jloCCil 


(34) 


where I is evaluated by summing up all local sources within 
the cell: 

jlccal ^ r° S{T')exp{-T')dT'. (35) 

Jo 

ALI is crucial for scattering-dominated problems whenever 
e < O.OI - see convergence tests in §3.2. Even with ALI, we 
often need to run hundreds of iterations before the solution 
settles down and converges. This is feasible with SC, which 
is quite fast, but not practical with the much slower LC. The 
convergence criterion that we use is: 

max ^ (36) 

where the convergence level Sc is preset to a desired accuracy 
level (typically Sc = 10“^). One caveat with this criterion is 
that it only indicates the level at which HERO is internally 
satisfied with the solution (i.e. that HERO converged onto a 
solution). HERO can still converge onto the wrong solution 
due to systematic errors caused by the numerical setup (i.e. 
the assumption of linearity for how S varies with r, or the 
assumption that /(fl) has linear dependence on across 
neighboring cells). We find that in practice, these systematic 
effects lead to the HERO solution being biased at the ~ 5% 
level when compared with exact analytical results. 
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2.5 Ray tracing 


2.1.2 2D Angles 


At the end of the day, we wish to calculate images and com¬ 
pute synthetic spectra from our radiating objects. This is 
accomplished by means of raytracing from some distant ob¬ 
servation plane. We consider a large number of rays arriving 
at the observer and trace each backwards in time to de¬ 
termine its trajectory and point of origin. Then, using the 
radiative transfer equation and the source function as calcu¬ 
lated via the LC method, we compute the observed intensity 
of each ray. Finally, we combine the rays to construct the 
observed image of the disc as well as the spectrum of the 
source. These synthetic results, which are easily calculated 
for different viewing angles, can be directly compared to ob¬ 
servations of accretion discs. 

Photon trajectories are handled by numerically evalu¬ 
ating the null geodesic equation Eq. 19 (a second order PDE 
in X and p): 


(fx^ 

dX^ 

with 


dx'^ 

d\ ■ 


(37) 


The emission seen at the observer plane is obtained by sum¬ 
ming up the emissivity along the ray paths (same exercise as 
Eq. 27). This yields a grid of ray intensities (’’image” of the 
system) as projected on the observation plane. Integrating 
all the rays generates the final observed spectrum via 

F^ = J hcosOdn, (38) 

where ^ ~ 0 is the normal incident angle of the ray at the 
observer plane, Ijy is the intensity at the observer plane, and 
dO, is the angular size of a single pixel at the observer plane. 


2.6 Frequency Discretization 

HERO is set up to handle frequency dependent opacities. 
Both the computational cost and memory requirement scale 
linearly with np the number of frequency bins. All radiative 
quantities (opacities, intensities, source functions) are calcu¬ 
lated for the same discrete set of frequencies. The code could 
be set up to handle group mean opacities and emissivities 
with their appropriate quadrature weights to handle more 
complex line emission problems. In problems where there is 
frequency coupling (i.e. gravitational/doppler shifting) the 
redistribution of photons is handled by linear interpolation. 
Our typical choice for the frequency grid is a total range 
of 6 decades with a resolution of 10 points per decade (60 
frequency bins). 

2.7 Angular Discretization 

The angular setup of our code differs based on the dimen¬ 
sionality of the problem that is being solved. 

2.7.1 ID Angles 

In ID, we solve plane parallel slab problems which are ac¬ 
tually 3D problems that can be collapsed down to ID by 
invoking translational and azimuthal symmetry. We subdi¬ 
vide the 47r steradian sphere into equal solid angle slices in 0; 
thus, we set our angle spacing dO such that sin^d^ = 2/Na- 


In all the 2D test problems described in this paper, the rays 
are considered to live in a flat 2D space, so we do not consider 
the 3D solid angle at all. We choose our angle grid so as to 
cover uniformly the full 27t of the 2D plane. We use simple 
linear interpolation on 2D angles to handle the mixing of 
angles when rays travel from one cell to another. 


2.7.3 3D Angles 

The goal is to subdivide the Air sphere of solid angle into 
Na solid angle wedges. A few common approaches include 
bisecting octants, following a ’’spiral” that winds around 
the sphere from the two poles, or a special grid of lati¬ 
tudes/longitudes. All these methods approximately achieve 
the required goal, but they suffer from undesirable patterns 
of symmetry lines or ’’seams” on the surface of the sphere. 
To sidestep this issue, we obtain our angular grid via a nu¬ 
merical approach. 

The strategy that we employ is to deposit Na ficti¬ 
tious charged particles constrained to move on the surface 
of a sphere, with initial positions set by simple heuristics 
(equal proper length spacings in ^,0 coordinates). Then, 
we allow the particles to evolve over time under their mu¬ 
tual inter-particle electrostatic forces until they settle on an 
equilibrium configuration. This naturally produces a set of 
positions that are nearly equidistant from one another. We 
use these positions to specify the angular grid (see Eig 3 for 
an example with Na = 80). 

Eor relativistic problems where light bending introduces 
mixing between different angle bins, or for curvilinear coor¬ 
dinate grid setups that invoke angle grids fixed along the 
locally rotated unit direction vectors, we handle angular in¬ 
terpolation via linear combinations. Eor a ray with a given 
unit direction vector '0, there will always by a set of 3 an¬ 
gular points ( 01 , 02 , 03 ) on the Air sphere which defines a 
triangle enclosing 0. We find the angles (0i,02,03) by sim¬ 
ply locating the 3 angles with the largest dot product to 0 
in the angle grid. 

Once the three vertices of the triangle have been iden¬ 
tified, we decompose 0 as a linear combination of 0i, 02 , 
03 : 


0 = Ci0i + C202 + C303. (39) 

The coefficients ci, C 2 , C 3 (after renormalizing to satisfy ci + 
C 2 + C 3 = I) are then used as the interpolation weights. In 
rare cases (i.e. when 0 is slightly outside the boundary of 
the triangle formed by 0i, 02 , 03 ), the interpolation weights 
may be slightly negative. Eor stability purposes, we clip any 
negative weights to zero. 

In Eig. 3, we show how our linear interpolation scheme 
handles an example intensity pattern. Notice that in the 
slowly-varying equatorial regions, the interpolated result 
provides a much better match to the exact intensity distribu¬ 
tion than simply using a discretized version of the intensity 
field. The linear reconstruction does poorly when the radia¬ 
tion field fluctuates on angular scales smaller than a single 
beamwidth of the Na = 80 grid (e.g. see polar regions of 
Fig. 3). 


© 0000 RAS, MNRAS 000, 000-000 



8 V. Zhu, R. Narayan., Sgdowski A., & Psaltis D. 





Figure 3. A plot of how our interpolation scheme performs 
compared to the exact solution. Top: Exact / = cos(4^) cos(40) 
beam pattern, Middle: discretized version of top panel using our 
Na = 80 predefined grid, Bottom: / generated from linear inter¬ 
polation of the Na = 80 grid. 


In all the following examples, we run HERO decoupled 
from any hydrodynamics. HERO is a postprocessor to com¬ 
pute the radiation field lu and, if required, a self-consistent 
temperature Tgas-) given the optical properties (opacities, 
emission, absorption) of the medium. Other properties of 
the background fluid, specifically density and velocity, are 
assumed to be given and fixed. It is also assumed that the 
system is time-independent, so we effectively solve for the 
steady state solution for the radiation. 

We are particularly interested in problems where scat¬ 
tering dominates the opacity (as is often true with relativis¬ 
tic accretion flows), leading to diluted-blackbody radiation 
fields. Due to the complicated nature of scattering, many it¬ 
erations are needed for the radiative solver to converge. This 
makes scattering problems ideally suited to test the conver¬ 
gence properties of the code. We describe below a series of 
ID, 2D and 3D tests. 


3.1 ID Plane-parallel Grey Atmosphere 


Scattering problems are in general very difficult to solve due 
to the integro-differential nature of the scattering kernel in 
the radiative transfer equation (see righthand side of Eq. 1). 
As such, there are very few examples that have closed form 
analytic solutions. 

The ID plane-parallel isothermal slab with grey opac¬ 
ity turns out to be one particularly simple case with a well 
known closed form solution in the two-stream limit (c.f. §1 
Rybicki & Lightman 1979). There is also a semi analytic 
solution in the limit of infinite angles. Eor this problem, all 
quantities, T, /^, a, are uniform within the medium, leading 
to a solution that is only a function of the optical depth: 
dr — [n a)dz. Under the two-stream approximation, and 
using the Eddington approximation to close the radiative 
moment equations, the analytic solution for the mean inten¬ 
sity is given by 


J(t) = B 


exp(—v/Ser) 

1 + v^ 


(40) 


where B is the thermal blackbody source function, r is the 
total optical depth from the surface, and e is the photon 
interaction destruction probability given by Eq.7. 

Note the exponential transition in r that separates a 
thermal interior from a dilute-blackbody surface layer. In 
Eigure 4, we show solutions computed by HERO compared 
to Eq. 40. The r = 0 surface outer boundary condition was 
set to have zero ingoing flux, and the innermost r = 10^ 
boundary was set to J = H for some fixed B. The agree¬ 
ment between the numerical and analytical solutions is quite 
good (~ 1% errors) with a systematic trend of larger errors 
in low e systems. These low e atmospheres are highly scat¬ 
tering dominated and the induced systematic error is simply 
a consequence of the nonlocal nature of scattered radiation. 


3 NUMERICAL TESTS 

We wish to verify the following properties of HERO: 1) that 
it correctly solves for the radiation field using both short and 
long characteristics, 2) that it calculates a self consistent gas 
temperature taking into account the condition of radiative 
equilibrium, and 3) that it correctly treats the curved space- 
time of GR. 


3.2 Convergence Tests 

The previous plane-parallel atmosphere calculations for the 
high scattering cases (e < 10“^) required ALI to converge. In 
Eigure 5, We show the convergence properties of the code 
for a few example problems to emphasize the importance 
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Figure 4. J computed from HERO (dashed lines) compared with 
the analytic solution (solid lines) for a plane parallel scattering 
atmosphere described under the 2-stream approximation. Four 
values of e (see Eq. 7) are shown. 


of acceleration in obtaining the correct solution within a 
reasonable number of iterations. 

In the absence of acceleration, radiation information can 
only propagate a distance (5 t ~ 1 per iteration. Thus, the 
number of iterations needed for the radiation to pass from 
one end to the other is ~ Ttot iterations (corresponding to 
random walk diffusion with stepsize 6t = 1). 

The power of ALI rests in the fact that it boosts the 
single iteration range of influence from 5t — 1 to 5 t — Srceio 
To see why this occurs, consider a tiny perturbation of the 
source function AS for some cell. The effect of AS on a 
neighbouring cell is attenuated by an exponential optical 
depth factor exp(— At), where At represents the inter-cell 
optical depth. The ALI boost factor (c.f. Eq. 32 when e ^ 1) 
in the scattering dominated limit is ASagi = AAS'exp(AT), 
which exactly compensates for the exp (—At) attenuation. 
This allows the AS information to propagate to the next 
cell in just a single iteration, hence requiring Nceii iterations 
for information to traverse the medium from end to end. For 
most scattering problems, this leads to a huge speedup (since 
Ncells < rlt). 


3.3 Multiray Temperature Solution 

Here we describe a test of HERO’s ability to calculate 
equilibrium temperatures. In Figure 6, we compute a 
constant flux nonscattering atmosphere and compare the 
numerical temperature profile with the analytic result 
(from a table of the Hopf function, see Chandrasekhar 
1950). With sufficient angular resolution {Na > 16), the 
temperature profile computed with HERO is correct to 
within 1%. 



Figure 5. Convergence rates for highly scattering plane parallel 
atmospheres. Note that the standard Lambda Iteration (LI) pro¬ 
cedure saturates at a fixed error threshold for small values of e 
(when e 10“^). Accelerated Lambda Iteration (ALI) using Eq. 
32 converges for all values of e, even extremely small values. 

3.4 Test of Spectral Hardening 

As a test of HERO’s multifrequency capabilities, we now ex¬ 
amine a frequency dependent ID atmosphere. This problem 
is simple enough to admit an exact analytic solution for the 
spectrum at all depths within the atmosphere. Appendix B 
derives the depth dependent spectrum. 

We are particularly interested in examining how a 
strongly scattering atmosphere modifies the thermal emis¬ 
sion from deep within a plane parallel atmosphere. This phe¬ 
nomenon is known as spectral hardening, and acts to shift 
the apparent colour of the photosphere emission. The effect 
is particularly important in the context of modelling black 
hole accretion disc spectra (Shimura & Takahara 1995; Davis 
& Hubeny 2006), where the high degree of scattering leads to 
a shift in the thermal emission peak by a factor of ~ 1.5 — 2 
towards higher energies. 

Figure 7 demonstrates an extreme version of this spec¬ 
tral hardening effect for an e = 10~® plane-parallel medium. 
HERO produces spectra that agree to within 5% of the an¬ 
alytic solution for all energies and optical depths. Note that 
above the photosphere (t < Teff = l/\/3e 10^), the ra¬ 

diation spectrum differs quite substantially from the local 
thermal emissivity/temperature (dotted curves). 

For simplicity, the example shown here is with the two- 
stream approximation for grey (i.e. frequency independent) 
homogeneous opacities. We impose a constant flux boundary 
condition at t = 10^ and solve numerically for the equilib¬ 
rium temperature and frequency-dependent radiation inten¬ 
sity as a function of depth. There is excellent agreement 
between the code output and the analytic temperature and 
radiation profiles (Eq. B13 in Appendix B). 

3.5 Effect of a Heating Source 

As a final accretion disc-like ID test that utilizes most of 
the features of the code, we discuss a problem that includes 
a heating source. We consider a slab with a total optical 
depth of 2 X 10^ between two free surfaces. We assume grey 
opacities with e = 10“^ (a strongly scattering-dominated 
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Figure 6. Comparison of analytic solution (solid line) to numer¬ 
ical results (dashed lines) for e = 1 constant flux atmospheres 
modeled with different numbers of angles. The bottom panel is a 
zoomed in version of the top panel. The temperature solution in 
the interior (r ^ 1) matches perfectly to the analytic solution. 
The agreement at the surface improves as the number of rays is 
increased. 


disc). The material in the slab is heated at a steady rate 
of 10“^ (arbitrary units) per unit optical depth, so that the 
steady state flux from each surface is unity. We use 20 angles 
and 61 frequencies. 

We start the calculation with some arbitrary temper¬ 
ature profile [B — 10^ at all r in this particular example) 
and we run HERO until it converges. After 3000 iterations 
(this is many more than needed, but we wished to converge 
as much as possible for this test), the radiative transfer 
equation has a maximum error of 10~® and the radiative 
equilibrium equation a maximum error of 10“^. The results 
are shown in Fig. 8. As expected, the temperature profile is 
quadratic in r and symmetric with respect to the mid-plane 
(upper panel) and the radiation flux is linear in r and an¬ 
tisymmetric (middle panel). The lower panel shows spectra 
at various depths. In the deep interior of the disc, the ra¬ 
diation is nearly blackbody, but as we approach the surface 
(Teff = T\/3e < 1) there are signatures of spectral harden¬ 
ing. The radiation that escapes from the surface is distinctly 
hardened. 



Figure 7. Comparison of depth dependent spectra calculated us¬ 
ing HERO (solid lines) and the analytic solution (dashed lines) 
for a plane parallel thermal atmosphere. Note: this is a highly 
scattering test problem with e = 10“®, which results in signif¬ 
icant spectral hardening at the r = 0 surface (compare dotted 
blue blackbody curve with solid blue line). At all depths and fre¬ 
quencies, HERO is able to calculate the spectrum to within 5% 
of the true solution. The spatial resolution chosen for this calcu¬ 
lation was 25 points per decade in r, and 20 points per decade in 
u. The two-stream approximation was used for handling angles 
(Aa =2). 


3.6 2D Solutions and Ray Defects 

A well known limitation of radiative solvers with discretized 
angular zones is the appearance of “ray defects”. These are 
sharp linear features that arise because interpolation on the 
angular grid is imperfect. The presence of ray defects is an 
important motivator for us to use a hybrid scheme. Our long 
characteristics method operates independently of an angular 
grid since it adaptively chooses an angular grid and hence 
does not suffer from ray-defects. In the following sections, 
we compare our two solvers, SC and LC, for a few test 2D 
problems. Appendix A provides a more extensive discussion 
of the ray defects. 


3.6.1 Opaque Wall Test 


Ray-defects are particularly severe when there is a compact 
source of radiation (see Fig. AI in Appendix A). But there 
are noticeable effects even in the case of smooth extended 
sources. As an illustration, we analyze a simple 2D example 
where it is easy to compute the exact solution. In the top 
panel of Figure 9, we show the radiation field for an empty 
(zero opacity) box illuminated by a hot spot in the centre 
of the upper wall. The boundary condition at this wall is an 
opaque (r = oo) isotropic source function with a gaussian 
profile: 


S{Xc) 


= exp 


(X - Xc)^ 

w 


(41) 


where Xc corresponds to the centre of the box, and the width 
w is equal to 1/5 of the box size. Since the interior of the box 
has no opacity, the radiation field at any point can be eas- 
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Figure 8. Top: Self-consistent solution with HERO for the tem¬ 
perature in the interior of a uniformly heated ID slab. Middle: 
Numerical solution for the frequency integrated radiative flux H 
compared with the analytic solution. Bottom: Radiation spectrum 
Jjy compared with the blackbody spectrum at the local temper¬ 
ature Bjy for optical depths (from below) 0, 1, 10, 10^, 10^, 10^ 
from the surface. 


ily found by tracing rays backwards and finding the source 
function corresponding to the ray’s intersection point with 
the upper wall. We can thus calculate the mean intensity J 
at every point inside the box. 

In Figure 9, we compare the exact solution to the result 
from our two RT solvers. For the SC solver (middle panel), 
a clear wave pattern appears along each of the angle grid 
directions (we used Na = 20 in this test). These are the ray 
defects mentioned earlier. The severity of the defects grows 
as we approach the upper wall, because of the sharp break 
in the intensity distribution there. On the other hand, LC 
perfectly reproduces the radiation pattern, even in the far 
field limit (compare top and bottom panels of 9). 

3.6.2 Shadowing Test 

Shadowing is another classic test that is useful for spotting 
systematic biases in a radiative solver. Simple moment clo¬ 
sure schemes such as FED (Levermore & Pomraning 1981) 
or Ml (Levermore 1984) typically have problems resolving 
the correct shadow structure. For instance, FLD has trouble 
maintaining the coherency of shadows over long distances 
due to its diffusive nature, and Ml has issues dealing with 
multiple light sources. 

In Figure 10, we consider the shadow structure pro¬ 
duced by an optically thick square box illuminated from 
above by an isotropically emitting wall. The top panel shows 
the true solution, with a clear umbra and penumbra. The 
middle panel shows the result using HERO SC with a crude 
angular grid Na = 16. Note that the beam resolution of 
20^ is our typical angular resolution in 3D problems (80 
rays in 3D), so this is a realistic example of how SC would 
perform in 3D. Increasing the angular resolution makes the 
ray defects less pronounced - the top “exact” solution was 
produced with Na = 1000 using our SC solver. 

The bottom panel shows the result with the LC solver. 
Notice the high accuracy of LC in handling the radiation 
shadow pattern. This is a common theme in all of our tests. 
LC is always much superior to SC. 

3.7 3D Solutions 

We now shift our attention to 3D test problems, and begin 
by discussing ray-defects. All the ray defects discussed in 2D 
(previous subsection and also Appendix A) are present in 3D 
as well, especially in the case of cartesian grids. A spherical 
polar grid (the most natural choice for accretion problems) 
eliminates some problems by introducing ray mixing via the 
curvilinear nature of the coordinate system. However, this 
is at the expense of introducing a particularly serious defect 
for radial rays moving out from a central source. Specifically, 
if one applies the SC method blindly on a spherical grid, one 
will obtain a constant radiation energy density and flux at 
large radius instead of the inverse square law fall-off one 
expects. 

Dullemond & Turolla (2000) discuss a way to “fix” the 
radial beam problem. They slightly modify the radiative 
transfer equation along the radial direction such that the ex¬ 
pected inverse-square falloff is recovered. We build on their 
suggestion, except that, instead of modifying only one ray 
(the radial one), we treat all rays equally and apply an arti¬ 
ficial diffusion that, when coupled with a logarithmic radial 
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Figure 9. Comparison of radiation fields calculated using short 
and long characteristics with an exact analytic solution for an 
opaque gaussian wall emitting into vacuum (colour log J). From 
top to bottom: a) exact solution, b) SC solution {Nj^ = 20), c) 
LC solution. Note the appearance of systematic banding/waves 
in the SC case. These are ray defects. 



Figure 10. Comparison of shadow patterns calculated using 
short and long characteristics for an isotropically emitting top 
wall shining on a central opaque box. From top to bottom: a) 
“exact” solution = 1000), b) SC solution = 16), c) LC 
solution. Note the discrete levels (ray defects) in the SC shadow 
pattern. 
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grid, naturally produces an inverse square falloff in the flux. 
Details of our diffusion method are explained in Appendix 
§A2. The following tests as well as those in §3.8 employ this 
ray diffusion scheme. 


3.7.1 Ring Benchmark 


One particularly simple test problem is an axisymmetric 
opaque emitting ring in vacuum, for which it is straight¬ 
forward to calculate the radiation field at any point in space 
analytically. For an infinitesimally thin ring emitting at the 
equatorial plane, the total light reaching any position f is 
given by a 1-dimensional integral. 



(42) 


where C is a constant specifying the emission per unit length 
of the ring. In spherical coordinates, we have 


J(r, e,(j)) = J 


c 


J .2 j.f 2 _ 2rr' sin 0 sin 0' cos((/) — 0') 




(43) 

where r' = 0.5, 0' — cj)' G [0, 27r] are the coordinates of 
the ring in the setup described here. 

In Figure 11, we compare the results of SC and LC to 
the analytic result; the constant C has been appropriately 
normalized to account for the finite emitting area used in 
the LC/SC calculations. The LC calculation captures the 
radiation field perfectly {6J/J ^ 10“^), whereas the SC re¬ 
sult shows strong systematic ray-defect patterns. Note that 
to enforce a 1/r^ falloff of the radiation field, the SC calcula¬ 
tion employs our ray diffusion scheme (otherwise, the result 
would be significantly worse). 


3.7.2 Dusty Torus Benchmark 

Previously, in §3.1, we demonstrated that HERO correctly 
computes both the radiation field and the gas temperature in 
ID. Unfortunately, there are no analytic solutions available 
for nontrivial 3D problems. Therefore, we turn to a stan¬ 
dard benchmark problem that has been widely discussed 
in the literature and use this problem to numerically com¬ 
pare HERO with other radiative codes. The model in ques¬ 
tion consists of an axisymmetric dusty torus Pascucci et al. 
(2004) with density structure given by: 

p(r, z)= po- Mr) ■ MO (44) 

Mr) = {r/rd)~^ 

MO = exp{-Tr/4[z/h{r)f} 
h{r) = 2 d(r/rd)®/® 

The opacities are tabulated and correspond to 0.12jam sili¬ 
cate grains (Draine & Lee 1984). Scattering is assumed to be 
isotropic, dominating in the wavelength range 0.2 — 1.0/am. 
A stellar point source is located at the centre of the disc. 
This point source shines on the disc and dictates its ener¬ 
getics and temperature structure. In HERO, the radiation 
emanating from the central point source is set to the ex¬ 
act stellar solution at the innermost radial cells of the grid. 
The propagation of this radiation outwards is then handled 
by the short (or long) characteristics solver. To avoid being 
killed by severe ray-defects, we include the diffusive term in 


the short characteristics solver as described earlier (see §A2). 
We apply free outflowing radiation boundary conditions at 
the outer radius of the grid. 

Given the above boundary conditions, HERO solves for 
the radiation field everywhere, both inside and outside the 
disc, as well as the self consistent disc temperature, i.e. the 
temperature that satisfies Eq. 26. To benchmark the code, 
we consider the most difficult example presented in Pascucci 
et al. 2004, the case corresponding to r = 100. The upper 
panels of Eigure 12 compare the disc midplane temperatures 
computed by HERO with the benchmark models. Overall, 
the agreement is reasonable - the slight differences are likely 
due to ray defects. We recover the temperature structure 
to within 10%, with the worst cells being located in the 
optically thick disc midplane. 

The equilibrium temperature profile is highly sensitive 
to the amount of scattered light (which dominates over the 
direct stellar illumination near the disc surface by an order 
of magnitude). The good agreement between HERO and 
Pascucci et al. (2004) indicates that : 1) HERO correctly 
handles/redistributes the scattered light; 2) the temperature 
solver is robust. 

As part of this test, we also show in Eig. 12 how the 
short and long characteristics version of HERO perform in 
determining the self-consistent disc temperature. Short char¬ 
acteristics does a reasonable job everywhere within the op¬ 
tically thick parts of the disc, but systematically underesti¬ 
mates the radiation field and temperature in the optically 
thin regions. Generally, SG has difficulties propagating ra¬ 
diation towards the coordinate poles. Panel 2 of Eigure 12 
shows the resultant underestimated radiation temperatures 
(~ 10% error) and that the bias increases with increasing 0 
resolution. 

The long characteristics solver does not suffer from this 
systematic error and recovers the correct temperature pro¬ 
file within tens of iterations throughout the optically thin 
region. This confirms that it is generally a good idea to run 
a few iterations of long characteristics after the short charac¬ 
teristics solver has been run to convergence. We make a spe¬ 
cial point to emphasize the importance of the LG pass. One 
might be tempted to instead run a high resolution SG pass 
to pin down a more accurate radiative solution. However we 
find that high resolution does not reduce the systematic bi¬ 
ases that plague the SG method. Despite the computational 
expense incurred by the LG method, it is the only way to 
obtain an accurate solution to the radiation field. 

Einally, we also show in the lower two panels of Eig. 12 
integrated disc spectra. Agreement between HERO and the 
other benchmark codes is within 20%, which is comparable 
to the spread amongst the four independent codes discussed 
in Pascucci et al. (2004). The HERO spectra were computed 
by ray tracing from a distant observing plane, making use of 
the complete radiative solution (i.e. Su) obtained from our 
SG+LG hybrid solver. 


3.8 GR Solutions 

3.8.1 Light Bending 

HERO is designed to solve for the radiation field in general 
relativistic curved spacetimes. One important effect is light 
bending which is demonstrated in the test problem shown 
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Figure 11. Comparison of true solution (left) with the results obtained with SC (middle) and LC (right) for an axisymmetric emitting 
ring (colours represent log J). Note the effect of ray defects producing a “spider” pattern for SC. The LC calculation matches the exact 
analytic answer to within 0.1%. 


in Figure 13. J is computed with both the SC and LC in 
the 3 panels for a beamed light source located just outside 
the photon orbit (r = 3M) propagating in Schwarzschild 
space-time. We see that the regions with the highest in¬ 
tensity follow the expected curved trajectory. However, the 
beam is broadened substantially. This is a consequence of 
the hnite angular grid (number of angles Na = 80) used in 
the computation. 

The middle and right panels of Figure 13 show LC so¬ 
lutions to the same problem. The LC beam remains nar¬ 
row and coherent, agreeing very well with the expected be¬ 
haviour for free-streaming radiation at the photon orbit. The 
middle panel shows the LC result using an emitting source 
whose beam has been artihcially broadened to match the 
angular resolution of the SC Na — 80 angular grid. The 
right panel shows the true resolving power of LC, where the 
emitting source corresponds to a (5-function in angle. 


3.8.2 Disc Spectra 

To check that our handling of Doppler and gravitational 
redshifting is correct, we consider the problem of black hole 
accretion disc spectra. This problem has been tackled nu¬ 
merous times over the years using many independent codes 
(e.g. Li et al. 2005; Davis & Hubeny 2006; Kulkarni et al. 
2011; Zhu et al. 2012) and is a simple but useful benchmark 
test. 

We place an optically thick, geometrically thin disc 
around a Schwarzschild (a = 0) black hole radiating as per 
the idealized Novikov & Thorne (1973) disc model. In this 
problem, the disk emission is treated as isotropic thermal 
radiation with hux given by the NT model. For simplicity, 
we ignore any spectral hardening effects since we treat the 
emission in HERO as emanating from a single equatorial 
grid cell (i.e. for this test, we do not resolve the photon dif¬ 
fusion process that gives rise to spectral hardening). The 


inner edge of the disk is hxed at the ISCO nsco = 6 and 
the outer edge is located at r = 1000. 

We solve for the radiation held in the disc exterior using 
HERO in full GR. In HERO, the calculation of the disc 
spectrum is typically handled in two stages - 1) we hrst 
solve for the 3D radiation held above the disc using the 
short/long characteristics solver (this is not needed in the 
present example because we specify the disc emission prohle 
and assume vacuum outside the disc), and 2) we trace rays 
backwards (via Eq. 22) from a distant observer plane located 
to create a synthetic image of the disc. 

The HERO code solves for the full three-dimensional 
source function and radiation held within the 2M < r < 
lOOOM spatial domain of the grid. This information is then 
fed into a separate ray tracing subroutine. In the ray tracing 
stage, parallel rays are shot towards the disk distributed ac¬ 
cording to a squeezed logarithmic polar grid (the same setup 
as Kulkarni et al. 2011) from an observer plane located at 
r = 100, OOOM. The hnal spectrum is generated by integrat¬ 
ing the hux across the observer plane. 

The top panel of Eigure 14) shows the computed image 
for the particular example problem. Integrating over this 
image yields the observed disc spectrum, which we show in 
the lower panel. The disc spectrum computed with HERO 
agrees very well with a previous calculation by Kulkarni et 
al. (2011), lending conhdence that HERO correctly handles 
GR effects. 


3.8.3 Iron Line Spectra 

In a similar vein as for the previous test, which was based on 
the NT continuum spectrum of the disk, we also benchmark 
our code via a calculation of the emission prohle due to 
a monoenergetic line (e.g. Fe-Ka as seen in many Seyfert 
galaxies and some microquasars). Eor a system that emits 
only monoenergetic (5-function lines, the hnal integrated line 
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Figure 12. Comparison of HERO with the most difficult r = 100 dusty torus benchmark test of Pascucci et al. 2004. The top two 
panels show slices of the self-consistent temperature solution (top left: radial slices, top right: poloidal slices). We also show the SC 
solution evaluated for two different choices of grid resolution: Nq = (53,106). We find that the SC solution systematically underestimates 
the radiation/temperature field near the polar regions, with the effect becoming enhanced at high resolutions. This is a consequence of 
the angular diffusion scheme used in SC (see Appendix A2 for more discussion). The bottom two panels compare spectra at different 
inclination angles as computed from LC raytracing. 



Figure 13. Light bending test for our two radiative solvers using a narrow laser beam injected tangent to the r = 3 photon orbit. From 
left to right: a) SC, b) LC with an artificially broadened beam to match the SC angle grid size, c) LC pure. The solid black line shows 
the analytic result corresponding to the null geodesic. 
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Figure 14. Upper panel: Raytraced image of a razor thin disc as 
viewed by an observer located with inclination angle i = 60. Col¬ 
ors indicate log of the frequency integrated intensity. The asym¬ 
metry is due to the doppler effect combined with gravitational 
redshifting. The central arc feature corresponds to a secondary 
image of the accretion disk that arises from strong gravitational 
lensing about the black hole. Middle panel: Integrated disc spec¬ 
trum as computed by HERO (red points), compared to the result 
with the Kulkarni et al. 2011 code (blue line). Bottom panel: 
fractional errors between HERO and Kulkarni et al. 2011. 




Figure 15. A comparison of the Fe-Ko; line profile as computed 
by HERO and RELLINE (Dauser et al. 2013) for nonspinning 
(blue, a/M = 0) and nearly maximally spinning (green, a/M = 
0.998) black holes. Line profiles are computed for different viewing 
inclinations: top panel shows i = 30, whereas bottom panel shows 
i = 60. 


profile depends only on the geometrical properties of the 
system (i.e. redshifting from doppler/gravitational curvature 
and lensing effects from the Kerr metric). 

We consider a Keplerian accretion disk in the Kerr met¬ 
ric whose line emission is modulated by a power law emissiv- 
ity profile F{r) oc r~^. We also set the domain of the disk to 
span from rin — risco out to Tout = 400. Figure 15 shows 
the resultant spectra computed for two different choices of 
BH spin and two different choices of viewing angle. We find 
good agreement between HERO and the benchmark code 
RELLINE (Dauser et al. 2010, 2013), with the largest dis¬ 
crepancies occurring in the low energy red tail of the line 
profiles. 

The setup in HERO is identical to that of the previous 
NT disc tests, except that a monoenergetic 6.4 keV line was 
used as the local emissivity instead of a thermal continuum 
source. In addition, we use much higher resolution for the 
raytracing grid since the low energy red tail of the line profile 
depends sensitively on how well the inner edge of the disc is 
resolved. 
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Figure 16. A calculation of the returning radiation for a razor 
thin accretion disc around a moderately spinning a* =0.9 black 
hole. The standard Page & Thorne (1974) luminosity profile is 
used to set the outgoing flux/rays from the disk plane (blue). 
We compare the returning flux as calculated by the two radiative 
solvers in HERO (red, green) with the result from Cunningham 
(1976) obtained via relativistic transfer functions (black). The 
black dashed line is an extrapolation of the values from Cunning¬ 
ham’s data table. 

3 . 8.4 Returning Radiation 

Another classic accretion disk problem is that of comput¬ 
ing the returning radiation due to relativistic light bend¬ 
ing around the black hole. For this test, we setup a razor 
thin accretion disk that emits according to the standard thin 
disc luminosity prohle (Page & Thorne 1974) and measure 
the amount of returning radiation incident on the disk. The 
goal is to benchmark both our SC and LC radiative solvers 
against the solutions of Cunningham (1976), who tackled the 
same problem by means of relativistic transfer functions. 

In HERO, we model the returning radiation problem 
with a grey calculation (1 bolometric frequency bin) on an 
axisymmetric polar grid with (ur,n0)=(6O,3O) restricted to 
the upper half plane. The spatial grid was set as uniformly 
spaced in angle and log(r). Boundary conditions invoked are: 
reflecting for the polar axis (to account for light that passes 
through the pole), constant flux injection at the equatorial 
disk plane according to (Page & Thorne 1974), and zero 
incident radiation at the inner and outer radial boundaries. 

In hgure 16, we plot the incoming radiation flux at var¬ 
ious locations above the disk and compare to the solution 
of Cunningham (1976) for a moderately spinning a* = 0.9 
black hole. We And that SC systematically overestimates the 
amount of returning flux at large radii, presumably caused 
by the angular interpolation bleeding some of the outbound 
radiation into inbound rays. LC on the other does not expe¬ 
rience any strong systematic biases, but suffers from a lack 
of resolving power at large radii since at these large dis¬ 
tances, the LC ray grid can miss the inner photon ring that 
is responsible for most of the returning radiation. 

3.8.5 Vacuum Test 

As a hnal test of the general relativistic capabilities of 
HERO, we turn to the problem of light propagation in vac¬ 


uum from an opaque spherical shell. In Schwarzchild geom¬ 
etry, it is a simple exercise to solve for the apparent angular 
size of a constant radius shell as viewed by an observer at 
some other (larger) radius. Eurthermore, if the surface of the 
shell radiates isotropically like a blackbody with a constant 
surface temperature, one can calculate the radiation quanti¬ 
ties (i.e. J, E, L) by simply integrating the constant intensity 
across the apparent solid angle subtended by the shell. Eor 
instance, the radial prohle of luminosity has a particularly 
simple form, scaling directly with gravitational redshift as 


Llocal — Loo(l + 2 :)^ (45) 


where Liocai = Eiocai47rr^ is the total luminosity as mea¬ 
sured by an observer at radius r and Loo is the luminosity 
at inhnity. In Eigure 17, we show the results for the luminos¬ 
ity as calculated by our general-relativistic short and long 
characteristics solvers. We normalize the luminosity by the 
redshift factors so the analytic solution simply corresponds 
to a flat horizontal line (i.e., constant luminosity as mea¬ 
sured at inhnity). 

In general, we And that our short characteristics solver 
systematically underestimates the luminosity prohle at large 
distances. These tests were carried out with the diffusion 
prescription described earlier, therefore if the shell radiated 
purely in the radial direction, SC would by construction give 
the correct answer. Here the surface radiates isotropically, 
and there is a deviation from the true answer because of 
angle interpolation and diffusion. The various SC curves in 
Eigure 17 correspond to different choices for the radius of 
the inner emitting shell. The high degree of similarity for ah 
choices of inner radius (i.e. ranging from from highly rela¬ 
tivistic ro = 3 in units of GMjc? to nonrelativistic ro = 10^) 
implies that the SC bias is independent of relativistic effects. 
It arises purely from the discretization of the spatial and an¬ 
gular grid. While the bias is not negligible - the luminosity 
is reduced by a factor of 2.5 at large radius - note that the 
luminosity would be a factor ~ 10® too large if we did not 
include diffusion. As with the other test problems, LC gives 
an essentially perfect answer. 

Panel two of Eigure 17 shows another interesting phe¬ 
nomenon, viz., there is a strong dependence between the SC 
bias on the choice of spatial grid resolution. As we increase 
^-resolution (keeping the r-resolution Axed), the luminosity 
prohle of our isotropically emitting shell exhibits a stronger 
bias to lower values. This effect is ultimately caused by the 
angle diffusion scheme that we employ (c.f. Appendix A2). 
The diffusion is tuned such that the radial ceh-to-ceh atten¬ 
uation of light recovers the inverse-square law. If the cell as¬ 
pect ratio in r —^ is too rectangular, then the diffusion has a 
strong directional preference. Light rays traveling along the 
short dimension of the cell are overattenuated since they hit 
cell boundaries more often than rays traveling in the long- 
dimension (the degree of diffusion is directly proportional to 
how often the light ray traverses cells within a given spatial 
distance). Based on our tests, grid cells should ideally have 
a square aspect ratio to minimize the error, and any ratio 
in excess of 2:1 should be avoided. 
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Figure 17. Variation in radial luminosity profiles as a function 
of a) relativistic effects (leasing, redshifting) and b) grid-^ resolu¬ 
tion. Note that the luminosity plotted is the corrected luminosity 
as measured by an observer at infinity. The analytic solution is 
therefore a constant L/{\ + z)^ for all radii. 


4 SUMMARY 

We have described in this paper HERO, a new general rel¬ 
ativistic radiative transfer code. The primary aim of this 
code is to model the radiation field in accretion flows around 
black holes. The unique features of HERO are: 1) a hybrid 
short/long characteristics radiative solver that enables ac¬ 
curate and fast modelling of complex anisotropic radiation 
fields; 2) implementation in a general relativistic framework 
taking into account the effects of light bending, doppler 
beaming, and gravitational redshifting. 

HERO is written as a post-processing code decoupled 
from the hydrodynamic evolution of the fluid. It computes 
the time-independent radiation field assuming a given fixed 
background fluid structure. Strictly speaking, this approach 
is valid only for problems where the fluid velocities are small 
compared to the speed of light (i.e. nonrelativistic flows). 
Alternatively, and this is the primary application we have 
in mind, it could also be applied to time-steady relativistic 
problems. 

We provide a detailed explanation of the long/short 
characteristics method used to solve for the radiation held 
and our approach for solving the self-consistent gas tem¬ 


perature. To verify that HERO produces physically correct 
answers, we have performed a comprehensive set of tests de¬ 
signed to examine the code’s convergence properties, accu¬ 
racy, and capability to handle multidimensional relativistic 
problems. We conhrm the well known result that 2D and 3D 
problems with compact sources suffer from signihcant ray 
defects in the far held when analyzed with the short char¬ 
acteristics method. We present an approximate hx which 
mitigates the effects in the case of a 3D spherical grid. How¬ 
ever, for accurate results, it is necessary to switch to a long 
characteristics solver which is unaffected by ray-defects. 

As the subject of a follow-up paper, we intend to apply 
HERO to radiative MHD simulations of accretion discs - 
particularly simulations undergoing super-Eddington accre¬ 
tion, where radiation feedback strongly dominates the dy¬ 
namics of the how. Using HERO, we will investigate the 
integrated spectra to see the role that self-shadowing and 
irradiation plays in these systems. This application requires 
a Comptonization module which will be described in our 
next follow-up paper. 
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APPENDIX A: RAY DEFECTS 

The method of short characteristics constitutes the primary 
workhorse of our radiative solver and is used to generate 
a good first approximation to the radiation field. However, 
“ray-defects” are a well known limitation of SC, which is 
why we need to follow up SC with the more accurate LC 
method. Here we discuss some of the properties and explore 
the cause of ray defects. Figure A1 shows a few examples 
of ray defects in 2D arising from point source emitters. The 
defects manifest as unphysical beam-like patterns far from 
the emitting source. 

Point sources generate the strongest defect pattern so 
we use them in the following discussion to illustrate the main 
issues. Point-source-like emission does appear in the prob¬ 
lem of black-hole accretion discs, e.g. the hottest innermost 
region of the disc shines extremely brightly and due to its 
compact spatial scale, acts like a point source at large dis¬ 
tances from the centre. 

For a fixed angular grid, the point source ray defects 
form a series of radial beams that reflect the underlying 
structure of the angular grid (see Figures A1 for a few ex¬ 
amples in different coordinate systems). Beam collimation is 
enhanced for rays travelling in directions where neighbour¬ 
ing grid cells cover a smaller angular size. This ’’grid-lattice” 
effect is best seen in panel c) (sheared box) in Figure Al, 
where the thinnest beams are those travelling to the upper- 
right/lower-left sectors (i.e. the directions where the neigh¬ 
bour points as defined in Figure 1 span the smallest angu¬ 
lar extent). The same, but less pronounced result is seen in 
panel d) of Figure Al (polar coordinate system). Here, the 
rays pointed radially inward suffer less dispersion than their 
radially outward counterparts, again for the same reason as 
in the sheared box case (tighter angular packing occurs for 
the neighbour cells at smaller radius). 

Finally, ray defects are particularly enhanced for rays 
that directly intersect neighbouring cell centres (as an ex¬ 
ample, note that the 45° rays exhibit overwhelmingly strong 
defect patterns in our cartesian setup in the top two panels 
of Fig Al). 

Ideally, we would like to represent the discretized radi¬ 
ation pattern as a smooth field instead of a superposition 
of laser beams. After much trial and error, we have arrived 
at an angular diffusion based solution for spherical log polar 
grids. To motivate our final solution, we first examine the 
root cause of the defects. 


Al Mathematical Origin of Ray Defects 

Ray defects ultimately arise due to the local nature of the 
short characteristics solver. Methods that operate using just 





Figure Al. Comparison of the radiation solution for a point 
source computed by short characteristics for different grid choices 
{Na = 20). Colour indicates log( J). From top to bottom: a) carte¬ 
sian grid (100x100); b) same as a), but using quadratic interpo¬ 
lation; c) sheared cartesian (100x100); d) Polar grid (60x100). In 
all cases, we place a delta function source function that emits 
isotropically. The pencil beams that appear are the ray-defects - 
physically, the far-held radiation pattern should be spherical and 
smooth but the limitations of treating radiation only locally in 
the SC method results in the thin beams (see discussion in §A1). 


© 0000 RAS, MNRAS 000 , 000-000 



















HERO: A 3D General Relativistic Radiative Code 21 


local propagation of light develop a “Pascal’s Triangle” char¬ 
acteristic beam pattern with increasing distance. If the radi¬ 
ation from a source cell is split with some linear combination 
to its neighbouring cells and the same method is applied 
uniformly for ah cells, then the propagation of light has the 
following characteristic shape - consider the case where the 
mixing coefficients are [ 1 / 2 , 1 / 2 ]: 


1 



The hnal far held pattern simply corresponds to the 
weights of a discrete random walk. This shape is a Gaussian, 
with a characteristic width of ic oc y/d where d is the total 
propagation distance from the source. The beam far from 
the source will therefore have angular size w/d ex 1/y/d 0 

as d ^ oc. Thus, at a sufficiently large distance, the propa¬ 
gation of light via short characteristics will inevitably result 
in a series of laser beams. 

The defect problem can be further exacerbated if a ray 
happens to pass exactly through a neighbouring cell centre. 
In this case, the mixing/interpolation coefficients are [1, 0] 
and the propagation diagram reduces to: 


1 



1 0 0 


Notice that for this limiting case the beam pattern does 
not spread at ah, resulting in a zero width beam of constant 
intensity at all distances. This kind of defect is particularly 
disastrous in the case of spherical coordinates since it af¬ 
fects ah radial rays. Thus, any compact light source near the 
centre of the spherical grid will develop serious ray-defects 
at larger radii. The effect is quite devastating if there is a 
strong point source at the centre. The radial ray defect acts 
to force a constant intensity on all radial rays, resulting in 
constant radiation energy density independent of distance 
rather than the expected inverse square fahoff. The natural 
way to correct non-spreading beams is to introduce some 
degree of artihcial broadening (i.e. force the mixing coeffi¬ 
cients to have some minimum hoor value). We describe our 
approach in detail in the next section. 


A2 Ray Defect Correction Schemes 

The most obvious approach to mitigate the impact of ray 
defects is simply increasing the angular resolution. However 
this approach is not feasible for complex problems with large 
spatial domains due to the unfavorable computational cost 
scaling with resolution. 

One idea to combat ray defects is to use higher order 
interpolation schemes in treating the propagation of radia¬ 
tion. Davis et al. (2012) found that quadratic monotonic in¬ 
terpolation works well and suggest that it be used over stan¬ 
dard linear interpolation. Unfortunately, this does not ad¬ 
dress the fundamental problem of beam spreading - higher 
order schemes are actually counterproductive in that they 
produce even thinner beam patterns and amplify the signif¬ 
icance of ray defects (Compare the top two panels of Figure 
Al for linear vs. quadratic interpolation) . 

Ultimately, the problem (as illustrated in §A 1 ) has to 
do with the lack of spreading for a single beam. This moti¬ 
vates us to implement a diffusive scheme to bring back the 
necessary amount of spreading and to set the intensity val¬ 
ues in such a way that an inverse square fahoff of radiation 
density and hux is recovered. 

Our approach is the following: every time we read out 
intensity values from our angular grid (i.e. the interpolation 
scheme described in §2.7), we apply a floor to our interpola¬ 
tion coefficients so that there is always some degree of mixing 
between nearby angles. That is, we modify the coefficients 
Cl, C 2 , C 3 in Eq. 39 to new values c^, C 2 , C 3 given by: 


Ci j 3w 

where w is the diffusion coefficient that controls the amount 
of angle mixing. A more physical way to understand this 
diffusion coefficient is to convert it into an attenuation factor 
for nearby cells. For a monodirectional beam (i.e. 7 = 0 for 
all but a single angle bin), the beam will attenuate by a 
factor a after propagating across a single cell, i.e., 

7"+i = ar, (A2) 

where from Eq. Al, we hnd: 

a = (1 + rc)/(l + 3ic). (A3) 

This motivates a method for choosing the appropriate value 
of w for the case of a log-polar grid (where the radial spacing 
is uniform in logr). Since we desire an inverse square law 
fahoff of the radiation held, we simply solve for a in the 
equation = 0.01 where Ndec is the number of cells per 

radial decade. For example if Ndec — 20 points per decade 
(our canonical choice), we have have a = 0.794 and w = 
0.149. 

In Figure A2, we show radial prohles of F calculated 
with the SC for a point source central star. Note that, in 
the absence of diffusion, the effect of radial ray defects is 
to produce an unphysical solution where F is constant with 
radius. 

In Figure A3, we show the results from this diffusion 
scheme compared to the exact solution for an off-centre 
equatorial ring source (e.g. the same setup as described in 
§3.7.1) for different choices of the spatial grid resolution. We 
hnd that our diffusive scheme suffers from systematic over- 
attentuation in the polar regions for spatial grids that too 
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Figure A2. Radial profiles of F for various choices of diffusion 
coefficient as calculated by short characteristics code on a spher¬ 
ical polar grid (nr,U 0 ) = (60,100). Given these grid dimensions, 
the diffusion coefficient must be set to u; = 0.149 in order to 
reproduce an inverse square falloff for flux. 

rectangular. We recommend keeping the grid cell aspect ra¬ 
tio as close to square (1:1) as possible in order to minimize 
this effect. 

The scheme described here is a generalization of that 
proposed by Dullemond & Turolla 2000. They applied the 
same factor a derived above except that they did it solely 
for the radial ray. For a central point source, which pro¬ 
duces only radial rays, their method and ours produce the 
same correct inverse-law behaviour of flux. However, when 
the source is extended, e.g., the spherical shell test problem 
described in §3.8.5, their method causes the intensity to fall 
too steeply. Our method also suffers a similar systematic 
bias, but it is less severe (see Fig. 17). In the limit of an 
isotropic radiation field, e.g., the interior of a blackbody en¬ 
closure, their method would still cause intensity to decline 
with increasing radius. Our diffusive method, on the other 
hand, would give the correct result, viz., no change in inten¬ 
sity. These differences are relatively minor. The main feature 
of both approaches is that they recover the inverse square 
law at least approximately. However, a major limitation of 
our diffusive approach is that it only works for the case of 
uniform log-radial grids. For more complex grids (especially 
those with nonregular arbitrary/adaptive meshes), the only 
generalizable solution is to invoke a few iterations of LC to 
resolve the inverse square falloff from each emitting source, 
fn particular, if the illumination comes primarily from a very 
limited set of point sources, then the explicit LC handling 
of just these point sources shoudn’t significantly impact the 
overall runtime. 


APPENDIX B: ANALYTIC ID ATMOSPHERE 
SPECTRUM 

In this problem, we set up a constant flux ID atmosphere 
subject to fixed grey opacities {pii, — const) with an ab¬ 
sorption fraction set to e = (i.e. a highly scattering- 

dominated atmosphere). The local source function is a com¬ 
bination of thermal emission and reflected light (c.f. Eq 6). 

The analytic solution to this problem can be easily ob¬ 
tained by considering the moments of the radiative transfer 


equation. We define the first few moments of the intensity 
field as 


1 



-1 

1 


Kv — 9 / ^ ludg-) 

J 

(Bl) 

-1 


where g = cos(^) is the angle cosine with respect to the 

plane normal. Using these quantities allows us 
moments of the radiative transfer equation as: 

to write the 

dHi/ , . 

Bi/j 

dTi/ 

(B2) 

II 

(B3) 


Combining Eqs. B2 + B3 and invoking the Eddington ap¬ 
proximation {Kiy = Ji//3) yields 


(B4) 

which is simply a constant coefhcient second-order inhomo¬ 
geneous differential equation in J^. This allows us to con¬ 
struct an exact solution using standard methods. The solu¬ 
tions to the homogeneous system are simply 

01 (r) = exp(v^T), 

02(t) = exp(—\/^ t). (B5) 

The particular solution Jp that satisfies the inhomogeneous 
system is given by^ 

OO T 

Jp{t) = J dr' + 4>2 {t) J dr', (B6) 

T 0 

where g = 3eB is the inhomogeneous function, and W de¬ 
notes the Wronskian, defined by 

(B7) 

dr dr 

= 2^/ik (B8) 

Putting everything together, the solution takes the form 

J ^ Jp + Cl 01 + C202, (B9) 

where the undetermined constants are set by the boundary 
conditions of the problem. At the r ^ oo inner boundary 
we expect J ^ H, so we must eliminate the exponentially 
growing mode by setting ci = 0. To set the surface boundary 
condition, we make use of the two-stream approximation and 
evaluate 

= (BIO) 


^ The particular solution is constructed via the “variation of pa¬ 
rameters” method 
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Figure A3. Ray defect pattern from an isotropically emitting ring for 3 different choices of grid resolution. Here we use a spherical 
spatial grid with locally defined ray angles. Panels: a) (ng = 26) square aspect ratio; b) (ng = 52) 1:2 aspect ratio; c) (ng = 104) 1:4 
aspect ratio. Notice the systematic reduction of radiation near the polar regions when ng increases and also the characteristic “spider” 
pattern arising from our choice of spatial grid interpolation. 




and enforce a surface boundary condition that is consistent 
with the Eddington approximation: 


which sets 


HjO) 1 

m Vs’ 



(Bll) 


(B12) 


The final step is to determine the thermal source {B(r)) 
that is consistent with our radiation solution from Eq. B9. 
This can be calculated using our two radiative transfer mo¬ 
ment equations. Since we have a constant flux atmosphere, 
dHidr = 0, which implies J(t) = B{r) from Eq. B2. We 
combine this with integrating the pressure equation (Eq. 
B3) to yield the full solution 


J(t) = 3Hr + J(0) 


(B13) 
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